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Abstract. We give a wide class of Lie-Poisson systems for which explicit, Lie- 
Poisson integrators, preserving all Casimirs, can be constructed. The integrators are 
extremely simple. Examples are the rigid body, a moment truncation, and a new, 
fast algorithm for the sine-bracket truncation of the 2D Euler equations. 



§1. Introduction 

Hamiltonian systems are fundamental, and symplectic integrators (Si's) have 
been increasingly used to do useful extremely-long-time numerical integrations of 
them. Wisdom [17] has used fast SFs to integrate the solar system far more efH- 
ciently than with standard methods; there are numerous examples illustrating the 
superior preservation of phase-space structures and qualitative dynamics by Si's 
[9,10,11,14]. Many Hamiltonian systems are not in canonical form but are most 
naturally written as Lie-Poisson systems, which generally arise as reductions from 
canonical formulations in more variables. Examples are rigid bodies, fluid-particle- 
field systems (e.g. magnetohydrodynamics or the Vlasov-Poisson equations [7,8]) 
and general relativity. 

Integrating Lie-Poisson PDE's first requires a truncation of the underlying Lie 
algebra, the most promising approaches at this time being the sine bracket ([5], 
associated with 5u{N) and applying to systems whose configuration space is the 
set of symplectic maps of a 2n-dimensional torus, such as 2D incompressible fluids) 
and, for localized distributions, the moment truncation of Scovel and Weinstein [15]. 
Then, to integrate in time, there are general methods which provide a symplectic- 
leaf-preserving Poisson map [4,6]. They are not only implicit but require evaluating 
functions like "e'^'^s" via Taylor series; hence they can be very slow. Although they 
can be of any order [2], the time-step must be kept small so that the implicit 
equations can be solved quickly by iteration — so the implicitness is no advantage. 

If Lie-Poisson integrators are to be as practical as standard symplectic integrators 
they should be simple and fast. To this end we describe the widest general class of 
such systems for which explicit methods are available. (Examples of such methods 
were first found by Ruth [13] for canonical systems, and by ChanncU and Scovel [2] 
for Lie-Poisson systems.) Our class includes the sine-bracket truncation of the 2D 
Euler equations — the sine-Euler equations — for which the new method is not only 
explicit but 0{N / log N) times faster than the standard implicit method. 
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§2. POISSON SYSTEMS & INTEGRATORS 

A Lie-Poisson system has a phase space M = R" 9 x, a Lie-Poisson bracket 
{F,G} — ^^Jij^- where Jij = cJ^^Xk (with c^^ the structure constants of a Lie 
algebra) and a Hamiltonian H: M^M. The dynamics x = {x, H} = JVH preserve 
the Poisson bracket; a Poisson integrator is one whose time-step map x^x' {x) also 
preserves the Poisson bracket. Symplectic splitting methods apply when the Hamil- 
tonian is a sum of terms each of which can be explicitly integrated for example, 
in a canonic;al Hamiltonian system, H = T{p) + V{q), which leads to standard 
symplectic integrators [13,14]. For Lie-Poisson systems, our methods depend on 
the following observation. We first form the set of all abelian subalgebras of M. 

Observation. Let 



S = {(7fe C {1, .. .n} : Jy = yi,j€ak} 

Then the Lie-Poisson system, with Hamiltonian H{(Jk) (i-e. H depends only on Xi 
with i <E (7k) is linear with constant coefficients. 

So if these linear systems can be solved exactly, an explicit first-order Poisson 
integrator for a Hamiltonian with p terms H = J2k^i Hk{(Tk) is 

if = exp(AtXi) . . . exp(AfXp) (1) 

where At is the time step and Xk = JVHk; that is, just integrate each piece 
of the Hamiltonian in turn. A second order symmetric method ("leapfrog") is 
^(iAt)^-i(-iAt) = 

exp(^^i) ■ • • exp(^Xp_i) exp(AfXp) exp(^Xp_i) . . . exp(^Xi). 

Methods of any order can be constructed by composing several such steps [2,9,16,18]. 
Clearly an arbitrary linear term could also be included in any of the Hk- 

Now S certainly includes the singleton subsets so we can immediately integrate 
Hamiltonians of the form H = (xk)- We choose not to break with tradition 

and include the rigid body as an example. 

Example: The rigid body. Here m S R"^ is the angular momentum, H = 
i^m^ _|_ ^ _|_ ^j^g j^-g a,lgebra is so(3) so 

— ms 7722 
J = I 7713 —mi 
— 7712 ITT'l 

This J has a Casimir C = |mp (the total angular momentum) so H = H — C/2/i 
has the same dynamics; this leaves only two terms in H. With Uk = f^kij^ ~ j^) 
and Rk{6) rotation by 6 around the axis 77ife, the map corresponding to (1) is 

m' = Rs{Atw3)R2{Atu)2)in 

— a "standard map" of the rigid body. The symplectic leaves |m|^ = const are 
clearly preserved. (We believe that this simple method (extended to higher order) 
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would improve the rigid body simulations of Austin et al. [1] who used the (implicit, 
non-Poisson) midpoint rule.) The same technique applied to the planar pendulum 
in the form 

■ -X2 



H 



1 



J = 




Xl 





length 



{x\, X2 coordinates in the plane, Xz = angular velocity; Casimir x\ + x^ 
of pendulum) does in fact give the standard map (see, e.g. [4]). 

Even if H is of the form J2 Hk{<^k), these methods are only practical if the 
resulting linear systems can be integrated exactly — diagonalizing them numerically 
would destroy their efBciency. This has not been a problem in the cases we have done 
to date: often the linear system is already upper triangular (after a permutation) 
and can be solved by back-substitution. 

Example: A moment algebra. Here M is the algebra of 2nd and 4th order 
moments of a distribution f{q,p)- Coordinates are {q"p^) = / q"jr\f{q.p) dqdp 
which we collect in a vector x. It could arise as a truncation of the evolution 
of f{q,p) by Hamiltonian vector fields on the plane, such as the Liouville equation 
f + {f, h}qp = 0, the 2D Euler equations, or the ID Vlasov-Poisson equations [3,16]. 
We have 
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-4X2 


-2X3 





-8X5 


-6x6 


-4x7 


-2x8 
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\ ip") J 




-8X7 
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As pointed out in [3], if H were separable in p and q as in the Vlasov-Poisson 
equations, i.e. H = T{x3,xs) + V{xi,X4), then an explicit Poisson integrator is 
possible, just as in the canonical case. But from the above observation is it clear 
that 

H = Hi{xi,Xi) + H2{x2,xe) + H3{x3,xs) + i?4(x4, X5, xe, X7, xs) 
also admits an explicit Poisson integrator. 

Extra nonlinear terms — In the canonical case one often has a nonlinear term in H, 
not a function of q alone, which can nevertheless be integrated explicitly (and con- 
veniently, i.e. without using nonelementary functions). One such is H{qi,q2 +P2), 
which arises in the nonlinear Schrodinger equation and in the Zakharov equa- 
tions [9] and which gives rise to the well-known splitting method for the nonlinear 
Schrodinger equation. This phenomenon is less common in Lie-Poisson systems 
because of the more complicated evolution of those Xj not appearing in H. But it 
can happen, e.g. let H = H{xiXs) in the above algebra: 



Xl = Axi 

X2 = 

Xs = -AX3 



xi{t) = e^*xi(0) 

X2{t) = X2(0) 

X3{t) = e-^*X3(0) 



4 



R. I. MCLACHLAN 



is 

xr 
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where a = 0:3(0), h = a;i(0), and A = 4a;2(0)-ff'(a6) is constant. This time-dependent 
Hnear system can be solved: let L{t) diagonalize its matrix; then because L~^L 

is timc-indcpcndcnt and tridiagonal, the system is transformed by L to a time- 
independent, tridiagonal one, which can be solved by standard methods. 

It is not always advantageous to split H into the fewest possible number of pieces. 
Splitting even further can lead to less coupled systems to solve and does not affect 
the error. 

The local truncation error of such methods can be easily derived using the 
Campbell-Baker-Hausdorff formula (see [9]). For a method of order I, it is a sum of 
all iterated Poisson brackets of order Z + 1 of the i/fe's, with coefficients depending 
on the order in which the H^s are taken in (1). It may look as if the error "increases 
with the number of pieces in H" ; but this would mean that the number of ODE's 
is also increasing, which will affect the error of any method. In fact, numerical 
simulations show that these methods have roughly the same truncation errors as 
the implicit Poisson methods. For example, for the rigid body, the above method 
was 1.6x less accurate but lOOx faster than the implicit method. 

§3. The sine-Euler equations 

Here we consider the motion of an inviscid incompressible fluid governed by 
the 2D Euler equations. The field variable is the vorticity uj{x,y), which is 27r- 
periodic in x and y. The phase space of vorticities is the dual of the Lie algebra of 
Hamiltonian vector fields in the plane [8] . We have 



J — ojydx — ojxdy 
H = J ipLodxdy where V^V = ~^ 



w = J{co) 



m 

6u) 



(2) 
(3) 



This algebra has an infinite number of Casimirs, which we may take as 

Cn = J 

In Fourier space eqs. (2,3) become 

^mn = m X nWm+n 



LOr, 



E 
117^0 



117^0 

m X n 



-co- 



m-t-n^— n 
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where m x n = TOin2 — m2ni and for real uj, uj-n = w* . There is a finite-dimensional 
truncation of (2) [5], the sine bracket 

Jmn = ^ sin(em X njuim+n mod JV (4) 

where e = 2'k/N and all indices are henceforth reduced modulo A'' to the periodic 
lattice —M<m<M where TV = 2M+1, and we shall take N prime for convenience 
(the extensions to nonprime and even N are straightforward). This J has N — 1 
Casimirs which approximate (7„ for 2 < n < N. Eqs. (4) specify the structure 
constants of su(A'') in an appropriate basis, see [5,19]. 
The most natural truncated H is 



M 

"-2 S 



nx,n2=—M 

giving the sine-Euler equations, first proposed by Zeitlin [19] : 

M 



ti^m — 7 , j — [5 '^m+n'^-n (5) 



1 sin(em x n) 



ni,n2 = — M 
n^O 



e 



where as in (4) all indices arc taken modulo N . As a numerical approximation 
of (3), these equations are only O(e^) accurate. However, once the vorticity has 
rolled up into small scales, standard (finite difference or spectral) approximations 
are not very accurate either, and do not possess the Poisson structure or conserved 
quantities of (5). Perhaps (5) is best regarded as an extremely interesting model of 
the 2D Euler equations until its properties are better understood. 

The approach outlined above now gives an C'(A''^ log A''), explicit Poisson inte- 
grator of (5), preserving all — 1 Casimirs to within round-off error — which is 
faster than the 0{N'^) needed just to evaluate the right hand side of (5). We take 

(Tk = {nk : < n < iV} e E. 

To split the Hamiltonian we need a set of modes that generates the entire lattice, 
e.g. 

K = {(0, 1)} U {(l,m) :0<m<N} 

(this is why it is convenient to take N prime — otherwise multiples of K don't cover 
the entire lattice) and then 

AT-l 

kGK n=Q ' ' 

We now solve the linear ODE's ut — JVHi^. Of course ojm = for m £ (Tk- 
From (5), the other modes decouple into 2M sets of N equations, of which we only 
need to solve M sets; the others are their complex conjugates (see Figure 1). The 
variables in each set are a translation of a^, say by j: let Zm = i^j+mk, then 

M 

— / ^ ayiZyji—n 
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where 

sin(enj x k) 
= £|nk|2 '^"■^ 

These ODE's are circulant and hence diagonalized by the discrete Fourier transform: 
let z = Fz where F is the DFT; then 

z = Az where A = diag(Fa) 

so the equations can now be integrated exphcitly. 

Summary of algorithm. 

• for k s do 

• for j = 1^^, . . . translation of k do [may be done in parallel] 

• with Zm = tOj+mk, Set z' = F~^e^*^Fz 

• copy (z')* into W-(j+„k) 

• end do 

• end do 

The whole procedure requires 3M{N +1) = ^{N'^ - 1) DFT's of length A^. As 
there are FFT's available for sequences of prime length A'' [12], the whole algorithm 
is ^^(A^^ log A^). Figure 2 shows the relative energy error (H{t) — H{0))/H{Q) for 
10^ time steps with At = 0.05, Af = 7, = 1 and H{0) = 0.75. As is expected 
from an integrator which is a symplectic map on the symplectic leaves of the phase 
space, the energy error docs not grow with time. The errors in the Casimirs, e.g. 
C2 = X^WnW-n are due only to roundoff error, and grow by about 5 x 10~^^ per 
time step. 

This explicit method exists because the only coupled terms in the Hamiltonian 
belong to the sets a^; it is fast because of the special form of the sine bracket in 
this basis. Preliminary simulations indicate that the vorticity does not roll up into 
the high modes and that the evolution can be followed for arbitrarily long times. It 
will be interesting to see what the implications of this model are for the ergodicity 
and statistical steady state of the 2D Euler equations. 

Acknowledgements. I would like to thank Clint Scovel for a careful reading of the 
manuscript, and for encouraging to me to struggle, to seek, to find and not to yield 
to Poissonology. 
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Figure 1. Mode splitting in the sine-Euler equations. Here N = 7 and M = 3. 
O shows modes in one term iJk in the Hamiltonian (k = (1, 1)). ■ shows modes 
which are coupled together in the linear system ui = JVi/k (here j = (2,0)); □ 
shows modes whose values are the complex conjugate of the ■ modes. 



Figure 2. Relative energy error to i = 5000, At = 0.05, N = 7. 



